function fI = chebinterp(varargin)

%Interpolates from spectral grids (either fourier or cheb) to arbitrary
%grid.  Fourier grid is assumed to be x = -pi + 2*pi*(0:N-1)/N and cheb
%grid is assumed to be x = cos(pi*(2*(1:N)-1)/(2*N)).  Input format is
%(f,xI,kind,dim).  dim input is set to 1 if dim is not inputed. kind is a
%string which is either 'fourier' or 'cheb'

% input: (xmin, xmax)
% want: (1, -1)

% manipulate: -((inp-xmin)/(xmax-xmin)*2 - 1)
%            = (inp-xmin)/(xmin-xmax)*2 + 1

f    = varargin{1};
xI   = varargin{2}(:);
xmin = varargin{3};
xmax = varargin{4};
if length(varargin) == 5
    dim = varargin{5};
else
    dim = 1;
    f = f(:);
end

if ((xmin > xI(1)) || (xmax < xI(end)))
  error(sprintf('interpolation grid outside of data domain:  data = (%g, %g), whereas  grid = (%g, %g)', xmin, xmax, xI(1), xI(end)));
end



xI = (xI - xmin)./(xmin-xmax)*2 + 1;

xI = xI(:);


sizef = size(f);
N     = sizef(dim);
NI    = length(xI);
Ndim  = ndims(f);

tol = 1e-2;

N = N - 1;
m  = (0:N)'; m = m(:,ones(1,NI));
x = cos(pi*m/N);
xI = xI(:,ones(1,N+1)); xI = xI';  
I  = find(abs(x-xI)<tol); 

cm = ones(size(m)); cm(abs(m) == N) = 2;  cm(abs(m) == 0) = 2;
 
Cmn = -(-1).^(m).*sqrt(1-xI.^2).*sin(N*acos(xI))./(xI - cos(pi*m/N))./cm./N; %this will be badly behaved near cheby points

if numel(I) > 0
    Cmn(I) = 0;
    
    for k = 0:N
        if k == 0 || k == N
            p = 2;
        else 
            p = 1;
        end
        Cmn(I) = Cmn(I) + cos(k*pi*m(I)/N).*cos(k*acos(xI(I)))./cm(I)/p*2/N;
    end
end    

        
F = multiprod(Cmn,f,1,dim);

order = [2:dim,1,dim+1:Ndim];

fI = permute(F,order);

    
end

function C = multiprod(A,B,dimA,dimB)

%multidimensional array product
%C(a1,a2,...aN,b1,b2...bM) = A(a1,a2,...k,...aN)*B(b1,b2...k,...bM)

sizeA = size(A);
sizeB = size(B);
NdimA = sizeA(dimA); %number of elements in dimA
NdimB = sizeB(dimB); %number of elements in dimB

if NdimA ~= NdimB
    error('contraction dimensions must agree')
end

%permute dimA of A to last dim and dimB to first dim
orderA = 1:length(sizeA); orderA(dimA) = []; orderA = [orderA, dimA];
orderB = 1:length(sizeB); orderB(dimB) = []; orderB = [dimB,orderB];
AA = permute(A,orderA);
BB = permute(B,orderB);

%reshape into matricies
NA = numel(A)/NdimA;
NB = numel(B)/NdimB;
A2 = reshape(AA,NA,NdimA);
B2 = reshape(BB,NdimB,NB);

s1 = sizeA; s1(dimA) = [];
s2 = sizeB; s2(dimB) = [];
sizeC = [s1,s2];

C = reshape(A2*B2,sizeC);

end
